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ABSTRACT 


This study examines the binary source correlation technique for deter- 
mining vertical profiles of the refractive index structure parameter, Ch. 
Theoretical intensity scintillation covariance functions and power spectra 
for atmospheric layers of depth Ah at a mean altitude of h are derived. 
These functions are related to the photoelectron counting statistics and the 
spatial covariance of photoelectron counts for binary point sources. A 
linear detector array in the exit pupil of an optical system is examined, and 
the effects of the detection process uncertainty are explored. A lower bound 
for the expected error of an experimental determination of the spatioangu- 
lar covariance of counts is derived. This error is then minimized via a least 
squares analysis using the redundant information available in pairwise 
multiple correlations of a signal from the detector elements of a ten element 
linear array. The refractive index structure parameter profile is then 
derived and found to be the undetermined coefficients of the spatial covari- 


ance weighting functions in the least squares analysis. 
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I. INTRODUCTI 


The optical remote sensing of refractivity fluctuations induced by 
turbulent mixing is of the utmost importance to studies of coherent 
electromagnetic (EM) radiation propagation in the atmosphere. A quanti- 
tative characterization of atmospherically induced refractivity fluctuations 
as a function of altitude is required in the design and system analysis of 
adaptive optical systems. 

Adaptive optical systems are systems that attempt to correct a perturbed 
wavefront via mechanical or electrooptic methods. These methods include 
systems with mechanically deformable optical surfaces and phase 
conjugation techniques using nonlinear optical materials. However, the 
implementation of these techniques to correct for atmospheric degradation 
of optical signal quality depends on a thorough probe of atmospheric optical 
effects. 

Typically, three quantities are used to characterize the atmospheric 
propagation. These are the isoplanatic angle (89), the spatial coherence 
length of the atmosphere (rg) and vertical profiles of the refractive index 
structure parameter Cn*. The purpose of this thesis is to study theoreti- 
cally a single technique, binary source correlation, with the intent of find- 
ing a robust optical technique to profile the refractive index structure 
parameter. 

Several different approaches to the profiling problem have been 
employed with varying degrees of success. These include both active and 


passive means. Active probes of atmospheric structure include acoustic 


sounders (up to approximately 1 km) [Refs. 1,2] and pulsed Doppler radar 

(3-20 km) [Refs. 3,4]. Also, direct profile measurements may be made using 
microthermal probes mounted on a balloon [Ref. 5]. Passive means of find- 
ing the propagation path variations of the structure parameter include 
direct inversion of the amplitude scintillation covariance function [Refs. 
6,7], spatially filtered apertures [Refs. 8,9] and the crossed beam or binary 
source technique. 

The crossed beam method was originally proposed by Fisher and 
Krause [Ref. 10] and also by Wang, Clifford and Ochs [Ref. 11]. Essentially, 
the crossed beam technique proposed that two optical sources and two 
receivers be arranged such that the beams from the sources cross at a point 
in space. The cross correlation of the receiver outputs allow the determina- 
tion of the turbulence characteristics at the beam crossing point. The 
crossed beam technique, hereafter known as the binary source method, was 
implemented experimentally by several French teams [Refs. 12-15] using a 
binary star for the source. This method had advantages not present in the 
aforementioned passive techniques. 

The greatest advantages of the binary source technique are high spatial 
resolution and numerical stability of the algorithm. This is in sharp 
contrast with techniques involving the inversion of an integrated profile. 
The inversion of the scintillation covariance is often unstable due to noise 
[Ref. 16]; and, the spatial filtering technique has limited spatial resolution 
due to the width of the altitude weighting functions [Refs. 8,9]. Addition- 
ally, an artificial source will be available on the relay mirror experiment 
satellite and it is hoped that these techniques will be directly applicable to 
profile determinations made using time delayed correlations from this 


source. [Ref. 17] 


II. THEORETICAL BACKGROUND 


To fully appreciate the binary source method requires an understand- 
ing of the theory of EM propagation in turbulent media. The seminal work 
in this field was carried out in the Soviet Union in the 1940s and 1950s and 
this work is summarized in Tatarski [Ref. 18]. Modern treatments of the 
propagation problem in terms of the formalism of scattering theory have 
yielded a better understanding of the propagation problem; however, this 
formalism has proven impotent in furthering the predictive power of the 
classical theory [Ref. 19]. 

The classical theory as developed by Tatarski and also summarized by 
Clifford [Ref. 19], relies on a statistical description of the properties of the 


turbulent medium. A summary of these properties follows. 


A. TURBULENCE IN GENERAL 
Turbulence is a property of fluid flows and these flows are governed by 
the Navier-Stokes equations. The flow is usually characterized by a dimen- 


sionless parameter called the Reynolds number, Re, with 
Re = vL/u. (2.1) 


The quantities v,L and pu are the characteristic velocity, length scale and 
kinematic viscosity respectively. For small values of Re viscous dissipation 
dominates and the flow is laminar. As Re grows in magnitude, then a 
critical value is reached (approximately 1000 in the atmosphere) beyond 
which fluctuations in the velocity field are no longer damped. For values of 


Re greater than Re critical the flow rapidly becomes chaotic. 
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Following Lumely [Ref. 20], there are several distinguishing character- 
istics of turbulent flow. These are: 

- irregularity 

- large Reynolds numbers 

- diffusivity 

- dissipation 

- three dimensional vorticity fluctuations. 

The characteristics listed above are essential elements of any turbulent 
flow. For example, the diffusivity of velocity fluctuations in a turbulent flow 
and the rapid mixing inherent in this process causes an increase in the 
rates of momentum and heat transfer. 

Turbulence is also characterized by three dimensional fluctuations of 
the vorticity; and, no flow of less than three dimensions is truly considered 
turbulent. The vorticity fluctuations take place at multiple spatial scales 
with each scale characterized by a different value of Re. 

The energy in a turbulent flow is ultimately dissipated by viscous 
friction. This dissipation converts the macroscopic kinematic energy of the 
flow to thermal energy; and, this occurs at the smallest or inner scale of 
turbulence. Fluctuations with a scale dimension less than the inner scale 
(1-10 mm in the atmosphere) are damped by viscous dissipation. 

The largest scale in a turbulent flow is determined by the boundary 
conditions. This largest or outer scale of turbulent motion is where energy 
is pumped into the flow. The energy then cascades from the largest to the 
smallest scale sizes adiabatically and the energy is dissipated at the inner 
scale. The range of scales between lo and Lo, the inner and outer scales, is 


called the inertial subrange. 


The inertial subrange is so called because the time scale associated 
with the inertial transfer of energy from larger to smaller scales is much 
shorter than the time scale characterizing dissipation; therefore, the 
energy transfer is essentially adiabatic. However, the Reynolds number is 
dependent on the scale and decreases as the scale decreases. When a scale 
is reached such that Re is less than the critical value then laminar motion 
results at that scale size and smaller. Using this fact and the definition of 
the Reynolds number, the inner and outer scales may be related (within a 


constant) by dimensional analysis [Ref. 18]. The resulting relationship is 
(1p/Lo)*8 = Re (22) 


with Re the Reynolds number of the flow as a whole. 
The existence of a well defined inertial subrange provides a handle by 
which the statistical analysis of turbulence is carried out. This approach is 


pursued further in the following sections. 


B. STATISTICAL DESCRIPTION OF TURBULENCE 

The Navier-Stokes equations are ill-posed since there are more vari- 
ables than equations; and, unsolvable unless certain restrictive assump- 
tions are made. Unfortunately, these same restrictions exclude the 
turbulent regime of solutions and a statistical approach is required to 
analyze the flow. 

The full statistics would consist of a knowledge of the multidimensional 
probability distribution for the velocity field. However, the statistical analy- 
sis 1s complicated by the fact that turbulence is a nonstationary process. 
Nonstationarity implies that the mean and higher order moments of the 


velocity field are also stochastic functions of time. Currently, the only well 
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founded statistical analysis of the velocity field statistics is due to 
Kolmogorov and is valid only for scales in the inertial subrange of fully 
developed turbulence. (Ref. 18] 

Though the ordinary moments of the velocity field distribution are non- 
stationary, a statistical analysis based on a second order quantity called the 
structure tensor is carried out by Tatarski [Ref. 18]. The structure tensor is 


defined in general by: 


Dij(r) = (IviQnt+n—-vilrs)) (vyjnt+r-vj(c) , (2.3) 


with the angled brackets representing an ensemble average. This quantity 
is approximately stationary for scale sizes in the inertial subrange. The 
structure tensor may be further simplified by three assumptions and these 
are: 

- local homogeneity 

- local isotropy 

- incompressible turbulence. 

Homogeneity implies that the structure tensor depends only on the dis- 
placement r and not a particular location in space. Isotropy implies rota- 
tional invariance or radial symmetry in that only the magnitude of r and 
not the direction is important. The assumption of local homogeneity and 
isotropy is a weaker assumption that reflects the fact that the structure 
tensor depends only on scales very close to r. With the third assumption of 
incompressibility, that is the divergence of y equals zero [Ref. 21], homo- 
geneity and isotropy yield the following expression for the structure 


function: 


D,r(r) = ([v-(r,4+r)-v7(r,)}2) (2.4) 


where D,;, is the projection of the structure tensor in the direction parallel 
to the displacement r. This function describes the statistics of homo- 
geneous, isotropic and incompressible turbulence. [Ref. 19] 

For values of r in the inertial subrange of scales dimensional analysis 


of the structure function yields: 
p= = Cy2r2/3 (275) 


with C,2 the velocity structure parameter [Ref. 18]. The structure param- 
eter is a measure of the intensity of the turbulence and it is directly related 
to the refractive index structure constant mentioned in the introduction. 
1. Atmospheric Refractivity Fluctuations 

To relate the velocity structure function to extensive variables such 
as the temperature or refractivity then conservative passive additives must 
be considered [Ref. 18]. Conservative passive additives are those quantities 
that have no effect on the statistical analysis of the turbulent dynamics; 
however, the conservative passive additive is transported and mixed by the 
velocity fluctuations. An example of a conservative passive additive is 


potential temperature 0, given by: 


6 = T(po/p)-286 (2.6) 


with T the absolute temperature in Kelvin, p the atmospheric pressure and 
Po the pressure at sea level. This quantity is independent of altitude. 
Following Clifford [Ref. 19], the potential temperature fluctuation is 


related to the refractive index fluctuation by: 


An = —79x10-6(p/62)A@ Gap 


with p the atmospheric pressure in millibars. Since @ is a conservative 
passive additive then so is the refractivity n. And as shown by Tatarski 
[Ref. 18], conservative passive additive structure functions obey the same 
power law as the velocity statistics except near the viscous convective 
subrange which is near the inner scale. [Ref. 20] Therefore, the structure 


function for refractive index fluctuations is given by: 
Dr(r) = Cy2r2/3. (2.8) 


The refractive index structure function is used to find the power spectrum 
of refractive index fluctuations and since Clifford's derivation [Ref. 19] is 


clear it is summarized here. 


2. Power Spectrum of Refractivity Fluctuations 


The refractive index at a point r can be decomposed into a mean and 
a fluctuating part 
n(r) = (n(r)) + n,(z). (2.9) 


If nj is an analytic function then the spatial frequency spectrum can easily 
be determined by finding the spatial Fourier transform of n;(r). However, 
since nj 1s a stochastic function then the spatial frequency decomposition is 


carried out using the three dimensional Fourier-Stieltjes integral 


n,() = | expGKepdN@&) (2.10) 
Forming the covariance of the refractive index fluctuations, 


B@+rr,) = (2, @+z,)n,@,)) , (22a 


the power spectrum is found by using the Wiener-Khinchin theorem. 


The Wiener-Khinchin theorem asserts that the power spectrum 
and covariance form a Fourier transform pair. Using this and by invoking 
homogeneity and isotropy, Clifford [Ref. 19] expresses the power spectrum 


of refractive index fluctuations as: 


1 
On°K 





| drrB (r)sin(Kr) (2.12) 


.¢) 


@ (K) = 


This expression is simplified further by use of the following relation 


between the structure function and the covariance: 


1 
B (0) ~ Br) =F Dr) (2.13) 


Using this relationship and in the limit that the outer scale and inner scale 


go to infinity and zero respectively the spatial power spectrum becomes: 
®,(K) = .083Cp2K-11/3 (2.14) 


Actually, this expression is valid only for spatial wave numbers such that 
2nL,—! << K << 2nl,-! with the limiting process used to derive it being only 
an analytical convenience. 
With these expressions in hand the first order theory of optical 

propagation in the atmosphere is addressed in the next section. 
C. FIRST ORDER THEORY OF EM PROPAGATION IN THE 

ATMOSPHERE 

The first order theory of the propagation of an EM wave in the atmo- 
sphere is essentially equivalent to first order scattering theory in which the 


Born approximation allows an approximate solution to the wave equation. 


This approach assumes single scattering of a monochromatic wave 
incident on the atmosphere and it is valid for relatively weak turbulence. 
For strong turbulence multiple scattering effects cause many of the statisti- 
cal quantities of interest to saturate to a constant value; however, except for 
long propagation paths through strong turbulence (1-2 km in the late after- 
noon) the first order theory proves adequate for many purposes. 

The first order theory assumes that the atmosphere has unit magnetic 
permeability and zero conductivity and that the incident EM wave has a 
sinusoidal time dependence. With these assumptions both Tatarski and 
Clifford form the scalar wave equation for the propagation of a component 
of the electric field. Clifford [Ref. 19] proceeds to solve this equation by the 
method of smooth and small perturbations and this solution is summarized 
here. 

1. Solution of the Wave Equation 

Since the electric field is proportional to the magnetic field, and 
assuming that depolarization effects are negligible, then Clifford shows 
that it is necessary to analyze only one component of the full vector wave 


equation. The equation to be solved is: 


V2E + k2n2E = 0 (2.15) 


with E one of the components of the electric field, k the wave number of the 
radiation and n is the position dependent atmospheric refractive index. 
The solution to this equation is found by expanding E in a series of 


decreasing terms: 


E=E,+E,+Eo+... (2.16) 
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with the zero order term corresponding to the unperturbed field, the first 
order term the single scattering term and higher order terms correspond- 
ing to multiple scattering. Retaining the terms to first order and using 


equation 2.9 then equation 2.15 is resolved into the two equations: 
V2E, + k2E, = 0 (2.17) 
V2E) + k2E] + 2k2n) E, = 0 (2.18) 


This separation is accomplished by substituting E = Ej5+E, and equating 
terms of the same order to zero while neglecting terms of order nj? or 
higher. At this point it should be noted that each term Ey in the perturba- 
tive expansion is assumed to be of order nj™. Following Clifford and 
Tatarski [Ref. 18-19] the unperturbed field is assumed to be a unit ampli- 


tude plane wave propagating in the z direction and Ey) becomes: 
Eo = exp (ikz) . (2.19) 


Substituting this into the source term of equation 2.18 then this equation 


becomes the nonhomogeneous Helmholtz equation for the perturbed electric 
field Ey: 
V2E) + k2Ey = -—2k2n,e1kz (2.20) 


The formal solution of this equation is well known and it is the convolution 
of the source term of equation 2.20 with the Helmholtz equation Green's 


function: 
ikl r—r' | 
cn Cs 2 ikz’ 
E@) = jie |. 6 rsa [ 2k n(x Je ] (2.21) 
with the integration over the scattering volume. 
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To simplify the formal solution further the Fraunhofer approxi- 
mation is made in equation 2.21. The Fraunhofer approximation assumes 
that A is much less than the size of the scattering refractive index inhomo- 
geneity so that radiation is scattered in a small forward cone. In fact, if the 
varying refractive index is perceived as a weak phase screen then 
qualitatively the angle of scattering will be a maximum for inhomogeneities 
on the order of the inner scale. Therefore, from simple diffraction theory 
the maximum scattering angle should be on the order of A/I> with A the 
wavelength; and, for inner scales on the order of millimeters then this 
angle is on the order of 10-4 radians. Specifically, the Fraunhofer approx- 
imation assumes that the vertical distance from scatterer to receiver, 
|z—z'l, is much greater than the transverse displacement, |p—p'!, from 
the z axis. Making the Fraunhofer approximation means replacing the 
term |r-—r'l in the denominator of equation 2.21 by |z-z'l and expanding 
the term |r-—r'! in a binomial series retaining terms to second order in the 


phase. The resulting expression for the perturbed field E; 1s: 








2 ikz 5 ' ee 
Phe 3, { ik(p-0 ny 
E,@)= e jo rexp ot C= (222) 


vol 


The physical interpretation of this equation is that of a spherical Huygen’'s 
wavelet emitted at the scattering weak phase screen. The fringes produced 
by the interference of the Huygen's wavelet with the unscattered plane wave 
are interpreted as the amplitude and phase perturbations observed in 
atmospheric optical propagation. 

For the task at hand the quantities of interest are the fluctuations of 


the amplitude and the intensity of the optical signal about the unperturbed 


values. To find the amplitude fluctuations both Tatarski and Clifford use 
the Rytov approximation which assumes that the solution of the wave equa- 


tion is of the form 


E = Aes ; (2.23) 


with A an amplitude and s a complex phase. Letting E, = Ajexp(isg) then 
the ratio E/E, becomes: 


iW es , 

+5 =] xP fi(s-s J (2.24) 
Oo O 

And, taking the natural log of equation 2.24 splits the ratio into real and 


imaginary parts as follows: 
( 1) ( A, ) . 
log \ 1+ E =log \ 1+ A. +1 (S-S,) (2:25) 


Since Ej is assumed much less than E, and also that Aj << Ag, then the 
logarithm can be expanded quite accurately in a Taylor series. Carrying 


out this expansion and retaining terms to first order yields: 


al 


A, ‘ 
E =a ti (S-S ) (2.26) 


Following Clifford this implies that the amplitude ratio is approximately 
the real part of equation 2.22 normalized by the unperturbed field E,. This 


yields the following expression for the amplitude ratio: 


A é iil (Ge) 
1_k aur k(p—p )] ~1* 
wk oS a r cos eee ae ~| (paz) (2.27) 





Tatarski shows that inherent in this approximation is the assumption that 


the amplitude perturbation is small over the distance of one wavelength of 
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the radiation. This is called the assumption of smooth perturbations and 
this condition is virtually always satisfied for optical radiation. 
Conventionally, a quantity called the log amplitude is defined such 
that: 
x(r) = In(A/A,) = Aj/Ag (2.28) 


Since from the weak scattering viewpoint the log amplitude, x(r), is equal to 
the normalized amplitude fluctuation Aj/A) then it will be used in all 
further calculations. 

To proceed further the stochastic variable n;(r') in equation 2.27 is 


expanded in a two-dimensional Fourier-Stieltjes integral: 


n,@) = | do(K,z') exp GKep) (2:29) 


with the expansion in a plane perpendicular to the propagation direction. 
The random amplitude do now becomes a function of z’, the location along 
the propagation path. Substituting this expression into equation 2.27 and 
carrying out the indicated integrations, Clifford obtains the following 


expression for the log amplitude: 


| er 
x(r) = Je¥2(sfardoase sin Bee (2.30) 


Using this expression the second order statistics of the observed amplitude 
fluctuations and the power spectrum of these fluctuations is derived in the 
next section. 
2: ovariance and Power Spectr of Log Amplitude FI ations 
The spatial covariance of the log amplitude fluctuations is found by 


using equation 2.30 and the relation: 
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B (p,z) = (x(p,+p,z) x*(p 2) (2.31) 


Following Clifford and inserting equation 2.30 into equation 2.31 the spatial 


covariance function is written as: 


o {{ iK-(@,+9-K'-p] oP K(z-2') 
B(p,z) =k \] € | dzdz sin [| 


Z " 
sin [A S2?] (dok,z)d0(K,2"*) (2.32) 


with the superscript * denoting the complex conjugate. The ensemble aver- 
age of the random complex amplitudes is expressed in terms of the two 
dimensional spectral density of the refractive index fluctuations, 


F,(K’, z—z'), by using equation 1.5 in Tatarski [Ref. 18]: 


(do(K,z')do(K,z")*) = 8(K-K)F_(K’,2'-2")d’Kd"K (2.33) 


with d(K—K') a three-dimensional delta function. The two-dimensional 
spectrum of refractivity fluctuations is defined using the three-dimensional 


spectrum ®, of equation 2.14 and is written as: 


co 


F(Kz-z" = | cko,a0.K, cosLK (z'-z")] (2.34) 


—oo 


Substituting equation 2.33 into equation 2.32, the observed spatial covari- 


ance of the log amplitude fluctuations becomes: 


200! 
Byo)= | axel {12 [| acac” sin [E22] 


sin [Re2")) F (K,z'-z") f (2.35) 
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This integral is the Fourier transform of the quantity in the curly brackets. 
Using the Wiener-Khinchin theorem the expression in the curly brackets is 
identified with the power spectrum of the log amplitude fluctuations. 

To proceed further the experimental geometry as it bears on the 
evaluation of the integrals for the covariance and power spectrum must be 
considered. The following sections evaluate the above expressions as they 


apply directly to the binary source method of structure constant profiling. 
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Ill. THE ANALYSIS OF THE EXPERIMENT 


The purpose of this calculation is to develop the binary source technique 
theoretically for application to a linear array of photosensitive or photo- 
emissive detectors in the exit pupil of an optical imaging system. The 
optical system is assumed to be well corrected with an approximately con- 
stant modulation transfer function over the range of spatial wave numbers 
included in the inertial subrange of turbulence. With this assumption, the 
intensity distribution in the exit pupil is related to the intensity distnbution 
in the aperture by a simple scaling. Therefore, this calculation examines 
the experimental geometry illustrated in Figure 1. The exact scaling rela- 
tions and corrections for a nonideal modulation transfer function will be 


developed in a later section. 


A. GEOMETRICAL CONSIDERATIONS 

The binary source technique determines the Cy? profile by finding the 
covariance of scintillations due to a finite number of atmospheric layers Ah. 
And, a consideration of the geometry of Figure 1 is used to simplify the 
analysis. 

The prototypical binary source of Figure 1 is a binary star with compon- 
ents labeled 1 and 2. Two portions of the wave fronts intercepted by the 
BP rector elements are illustrated and the following quantities are defined: 

d =1r1-r2 with d the detector element separation from center to center. 

Ad, the detector element width. 

h, the mean crossing altitude of the wavefronts as illustrated. 


Ah, width of crossing. 
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ha,hy; the maximum and minimum altitudes of crossing, respectively. 


8, the angular separation of the binary components 1 and 2. 


2. 1. BINARY 
SOURCE 


ALTITUDE 


Ah 


DETECTOR 





Figure 1. The Experimental Geometry 
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With these definitions and examining the geometry of Figure 1, two 


relationships between these variables are apparent: 


d= 6h (3.1) 
Ad = @Ah/2 (3.2) 


Equation 3.1 gives the connection between the detector element separation, 
d, and the crossing altitude, h. The second equation gives the relation 
between the detector element width and the width of the crossing at a mean 
height of h. Equation 3.2 is simply derived by considering the difference 
between the maximum and minimum crossing altitudes and using equa- 
tion 3.1. Note that equation 3.1 assumes that the detector array is parallel to 
a line segment joining the apparent positions of the binary source compo- 


nents. If this condition is not satisfied then equation 3.1 becomes: 
d = 60h sec 6 (ia) 


with @ the azimuthal angle between the array and the line joining the 
source components. Equation 3.1 is assumed valid for this calculation; 
since, the optical system can in principle be aligned and driven such that 


this relation holds. 


B. POWER SPECTRUM OF LOG AMPLITUDE FLUCTUATIONS IN Ah 
Given a picture of the experimental geometry the derivation of the 
power spectrum of log amplitude fluctuations for a layer Ah at an altitude of 
h is carried out using equation 2.35. This power spectrum is then used to 
find the expected spatial covariance of the log amplitude fluctuations and 
the related covariance of intensity scintillations via the Wiener-Khinchin 


theorem. 
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For a layer Ah, isotropy and equation 2.35 imply the following expres- 


sion for the power spectrum of the log amplitude fluctuations: 





Ah 
F (K.0)=k?| Jaz'az” sin KA 
ic ce 9k 
O 
Regea 
sin ’ Sa F (K,z = 7, ) (3.4) 


with Fy,(K,z—z') equal to F,(K,z'—z') the two-dimensional power spectrum 
of refractive index fluctuations [Ref. 19]. Also note that the origin of the 
coordinate system is at hg (z=0) and the integration is over the scattering 
volume or layer. 

To make full use of the symmetry in equation 3.4 the following variable 
change is made [Ref. 18]: 


= 7-7) meandas2ie= 2ebzp (3.5) 


Using the Jacobian of this transformation then, 


CVA Oz" 


V2 -1/2 
i] *t 0g og 
dzidat = |e | dédn = dédn (3.6) 
ge i a i 
on 7 an 


therefore, dz'dz' =d&dn. Using this variable change and the trigonometric 


identity sin(a+b)sin(a-b) = sin?a-sin2b with, 
a= Ken ; b=— (¢) (Gay 


then equation 3.4 becomes: 
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F (K,0) =k Jase, (K, §)L sin? eat zy -sin?(£2)] (3.8) 


The integrand is simplified further by noting Tatarski'’s observation on 
page 135 of Reference 18 that F,(K,&) rapidly approaches zero for K€ greater 
than or equal to unity. This derives from the fact that fluctuations in the 
refractive index separated by displacements greater than Lo, the outer 
scale, are uncorrelated; and, by the Wiener-Khinchin theorem the power 
spectrum vanishes as the corresponding covariance of fluctuations goes to 
zero. Also note that the assumption of smooth perturbations, A << lo, 
implies that k >> K over the significant part of the integrand of equation 3.8. 
Therefore, over the significant part of the range of integration then: 


Z 2 
Ke eK SL 
a << | and sin = ja (3.9) 
And, equation 3.8 simplifies to: 
2 2 (Kew 
F (K,0) =k {| dédn F _(K,¢) sin —) (3.10) 
én 


The evaluation of this integral depends on the assumption of a form for 
F,(K,&) and this form depends on the magnitude of Ah compared with the 
inner and outer scale. For the experimental geometry of Figure 1 the ten 
detector elements of the linear array will divide the total atmospheric 
volume into ten separate non-overlapping regions. The width of these 
regions will be on the order of kilometers for observations of the profile from 
the planetary boundary layer to about 20 kilometers. This is much greater 
than outer scales observed in the atmosphere [Ref. 2,19] and it is safely 


assumed that Ah >> Lo. 
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Under the conditions discussed in the previous paragraph the layer Ah 
includes many separate, uncorrelated regions of turbulence. Following 
Tatarski (Ref. 18] the power spectrum of refractive index fluctuations for 
sublayers of Ah on the order of the outer scale is still given by equation 2.14; 
however, the refractive index structure constant is now a function of n. 


Therefore, the refractive index fluctuation power spectrum takes the form: 


2 
F (Ken) = C_(n) f (KE) (3.11) 


and, it should be noted that: 
-11/3 () 
Jracenae =me0o3) K =! @ (CK) (3.12) 


with this expression completely analogous to equation 2.14. With an 
expression for the refractive index power spectrum available further 
progress in evaluating the integral of equation 3.10 is possible. 


Noting that equation 3.10 is even in € and symmetric in ny then substi- 


tuting equation 3.11 into equation 3.10 and employing symmetry yields: 


Ah/2 Qn 


Z 
2 2 24K (=n) ) 
F’ (K,0) = 2k | ancl sin (Ei | sacar 
8) 


oO 


Ah 2(Ah—n) 
2 2 K*(z—n) 
+ 2k | dC (n) sin (se f (K,§)d¢ (253) 
Ah/2 O 


for the log amplitude power spectrum. Over most of the region of integra- 


tion both 2n and 2(Ah-n) have magnitudes between the outer scale and Ah. 


As previously noted, the refractive index power spectrum rapidly 
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approaches zero for K& greater than or equal to one [Ref. 18], therefore the 
upper limits for the € integration are replaced with infinity to a high degree 
of accuracy. Replacing the upper limits of the € integrals of equation 3.13 
and carrying out the integrations using equation 3.12 yields the following 


expression for the power spectrum: 
Ah 


F (K,0) = ae fo (K)C* (n) sin” [Kew Jan (3.14) 
aE n n 2k 
8) 

Since the detector elements have a finite width then by equation 3.2 
spatial variations of C,* smaller than Ah cannot be resolved. In effect the 
detector observes an average structure constant for the layer. If the Cy? 
profile is known apriori then this average can be carried out in the Fourier 
transform or spatial covariance domain as an aperture average over the 
finite detector element area. However, since a measurement of the struc- 
ture parameter profile is the objective then it is the power spectrum that 
must be averaged over the layer depth. This average of equation 3.14 is 


indicated as: 


Ah 


(F (K,0)) = 2k"n Jos (K) (cn) sin” [xe an (3.15) 


0 


with the angled brackets denoting the spatial average over the layer. The n 


integration is trivially: 


(F 0) = 2k"n ©°(K)Ah (cq) sin® [en)) (3.16) 
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Past observations of C,? indicate that it fluctuates rapidly in space with a 


mean that decreases slowly with altitude [Refs. 2,19]. Since the spatial 


fluctuations of C,2 appear uncorrelated with the sin? weighting function 


then the average indicated in equation 3.16 is rewritten as: 
2 
2 2 . 2, K (zn) 
(F (K,0)) = 2k x @” (K)Ah (co?) (sin [AZ] (3.17) 


Carrying out the spatial average explicitly for the sin? weighting function 


yields: 
A 


h 
2 
1 |. 2[ K*@ 
(F (K,0)) = 2k"nAho® aK CRM) |i Lee dn (3.18) 
0 


or, for z, the observer coordinate, equals hg: 


(F (K,0)) = 2k2n aho® (K) (C2(m)) 
[1 (K-24) | (Kah) a 
mE a sin (> (3.19) 


with (Cy,2(n)) the average or effective value of the structure parameter for 
the layer at a mean altitude of h. This is the observed average power 


spectrum of the log amplitude fluctuations due to turbulent mixing in a 


layer Ah. 
The power spectrum used by the various French teams [Refs. 12-15] is 


in this calculation's notation written as: 


2 
; Kh 
F’ (K,0) = 2k? zsh? (K){C2(n)) sin? K*) (3.20) 
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Approximating the sin(K2Ah/2k) term in equation 3.19 by AhK2/2k (k >> K) 


then 3.19 becomes: 
2 
2 2 _ 2 K “j 
(F (K,0)) =) nAho?(K)(C2(n)) sin es (3.21) 


an expression identical to the French spectrum. The two spectra are illus- 
trated in Figure 2 and it is evident that the French spectrum is an excellent 
approximation to the spectrum of equation 3.19. Therefore, the spectrum of 
equation 3.21 will be used for further calculations since it is analytically 
simpler. 
C. THE SPATIOANGULAR CORRELATION OF INTENSITY 

FLUCTUATIONS 

The power spectral density of equation 3.21 can be Fourier transformed 
to yield the spatial covariance of the log amplitude fluctuations of a single 
point source due to a layer Ah. The solution of the wave equation in Section 
I] assumed a plane wave solution and as is well known from classical 
optics a point source at infinity produces plane waves at the detector. For 
the stellar sources of Figure 1, the transverse spatial coherence function is 
essentially unity for displacements of less than five meters as shown by 
Hanbury-Brown and Twiss for Sirius [Ref. 22]. Therefore, an assumption 
of initial spatial coherence of the source is excellent for spatial separations 
of less than five meters. However, the log amplitude covariance is not 
directly obtained by cross-correlation of the signal from the detector 
elements. 

For a real detector operating at optical frequencies the detector response 


is typically proportional to the intensity and it is the covariance of intensity 
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Figure 2. The French Spectrum (x), and the Spectrum of 
Equation 3.19 (-) 
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fluctuations that is probed by the detector system. The details of the 
detection process will be treated in later sections; however, the intensity 
scintillations in the aperture are independent of the detection process. The 
spatioangular covariance of intensity fluctuations due to the binary stellar 
source will now be derived and related to the log amplitude statistics. 
1. The Spatial Covariance Functions 

For calculational simplicity, consider the wavefronts and two 

detector elements of Figure 3. The normalized spatial covariance of the 


intensity is given by: 
I(r,) - UG@,) Moe) = Le) 


(cr,)) (Ie)? 


with the angled brackets denoting an ensemble average and I(rj;) the 
instantaneous intensity at rj. The sources labeled 1 and 2 are assumed 


independent; therefore, 
I(rj) = I1(rj) + Io(xj) G22) 


with the subscripts 1 and 2 denoting the sources. Note that conservation of 
energy requires that the long term average illumination in the exit pupil is 
uniform. Substituting this expression into equation 3.22 and rearranging 


terms yields: 
Ci(p) = ({1) + (,)) [Carp S1,(c,)) “- (81, (c,)81,(c,)) 


+ (oI) d1ye,)) + (31,0,)81,,)) (3.24) 
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Figure 3. The Wavefronts Intercepted by Detector Elements 
at rj, re Due to the Binary Source 
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with 51=I—(I). Defining the relative irradiance fluctuations as: 


roe (3.25) 





and substituting into equation 3.24 gives the following expression for the 


normalized spatial covariance of intensity scintillations: 


2 2 


a (1,) 
GG 8 (cow! | (fy) + (1 ya (1 Dy Foe ale Ge 1) 


ene 
(G+, 


The first two terms of equation 3.26 are the weighted autocovariances of 


C,(p) = 


> [4,(,) lye) + Ay) Lyfe)? (3.26) 


sources 1 and 2 with the weighting functions carrying the brightness or 
magnitude difference between the two stellar components. Since assuming 
homogeneity implies that the normalized autocovariances are equal then 


equation 3.26 is further simplified to: 


(1) + (Ly (1,) (I, 
SD) SEN (1c) 1) a ae ae) 
({1,) + (1,)) ((1,) + (1,)) 


[1,@,) 1,@,)) + 1,@,) 1()] (3.27) 


The first term of this equation is the normalized autocovariance of a single 
source and the second term is the cross covariances of source components 


1 and 2. 
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The intensity covariances in the observation plane are related to the 
log amplitude covariances by a relation due to Fried [Ref. 23], and since his 
derivation is clear and readily available it will not be repeated here. In this 


calculation's notation Fried's result is: 
Cy(p) = exp[4B,(p)] — 1 (3.28) 


Since the perturbative solution of the wave equation assumes weak turbu- 
lence and this implies that the log amplitude fluctuations are small then 


equation 3.28 is expanded in a Maclaurin series to first order. The result is: 
Ci(p) = 4Bx(p) (3.29) 


Both expressions for the spatial covariance of intensity are graphed 
in Figure 4 as a function of the log amplitude covariance. The first order 
approximation is selected for further calculations based on comparison 
with the experimental curve of Figure 4. As is evident from Figure 4 
experimental evidence indicates that the intensity covariance saturates to 
unity as the intensity of turbulence increases. Therefore, the first order 
approximation of equation 3.29 is actually a better approximation over the 
range of validity of the single scattering theory. [Ref. 24] 

Using equation 3.29 and examining the terms of equation 3.27 indi- 
vidually then the spatial covariance of intensity scintillations can be cast in 
a more transparent form in terms of the log amplitude covariance. For the 
first or autocovariance term of equation 3.27, direct substitution of equation 


3.29 yields: 


(1,@,)1,@)) = 4B, (p) = 4B, (p) (3.30) 


However, the cross covariance terms are not treated as straightforwardly. 
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Figure 4. The Normalized Variance of Intensity Scintillations Versus 
the Log Amplitude Variance for Various Forms for Cy], and 


Experimental Data Points (x) 
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The first cross term in equation 3.27 is expressed as a function of 


the log amplitude using equation 3.29: 


(1c, l(t) = 4B, ©) (3.31) 


This term is the cross covariance of the light propagating along the paths 
labeled by the numbers 2 and 3 in Figure 3. Using equations 2.31 to 2.33 and 
the assumption of statistical homogeneity then the cross covariance of log 


amplitude fluctuations is expressed as: 


BL av = | ex [iK-[(-z,) +p] F (K,0)d°K (3.32) 
The assumption of statistical homogeneity in the layer Ah implies that 
fluctuations in xj are mirrored by fluctuations in xg at the crossing. 
Further simplification of equation 3.32 is achieved by using equation 3.1. 
Substituting this relation between the mean altitude and the detector 
separation in equation 3.32, this cross term becomes the same function as 


the autocorrelation but with the origin shifted by @h: 


B. (p)=B (p-@6h);  r,-r, =— Oh (3.33) 
Xo 3 x 


An exactly analogous derivation for the second cross term yields: 


Oe Pid = B (p+@h) sae 1 6h (3.34) 


That is light propagating along paths 4 and 4 in Figure 3 have a maximum 
in the spatial covariance for refractive index inhomogeneities with a scale 


of 8h. 
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Substituting equations 3.33 and 3.34 in equation 3.27 then the spatial 


covariance of intensity scintillations is expressed as: 


(1,) +1), (1,) (I, 
) = + _ 4B (p) + —>—*— 


° 4[B (p-0h) + B (p+6h)] (3.35) 


iP 


with the log amplitude covariance given by the two dimensional Fourier 
transform of the power spectrum of equation 3.21. Also note that equation 
3.35 depends explicitly on 8 and the covariance is now called the spatio- 
angular covariance as a remainder of this dependence. 
2. Evaluation of the Spatioangular Covariances 
The normalized spatial covariance of equation 3.35 indicates that 
this intensity expression can be constructed from the log amplitude auto- 
covariance by appropriate scalings and corresponding shifts of the origin. 
Using the Wiener-Khinchin theorem the log amplitude autocovani- 


ance in the aperture plane is given by: 


Co 


B (p) = Jexplikep] F (K,0)d°K (3.36) 


—09 


with the power spectrum given by equation 3.21, and letting the average 
spectrum condition be understood. Note that (C,2(n)) is now identified as 
the mean or effective structure parameter in the layer at a mean altitude of 
h. Equation 3.36 is the two-dimensional Fourier transform of the power 
spectrum. However, the assumption of isotropy and an assumption of a 
circular aperture imply that equation 3.36 can be recast as a Hankel 


transform: 
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circular aperture imply that equation 3.36 can be recast as a Hankel 


transform: 


oo 


B (p) = 2n [F (KOM (Kp)KdK (3.37) 


0 
with J,( ) the zeroth order Bessel function. Substituting from equation 3.21, 


the autocovariance becomes: 


[oe] 


2 
B(p)=M|J (Kp) K °° sin°(32)aK (3.38) 


0 
with M=4n2(.033)k2Ah(Cy2). Note that this integral is valid only for the 
inertial subrange of spatial scales since the power spectrum is divergent at 
the origin. Therefore, the limits of integration are changed to K'=2n/Ly and 
k"=2n/l,, with 1, and L, the inner and outer scales respectively. 


In the limit that p goes to zero then equation 3.38 becomes: 


K" 


K*h 
B (0) =M fi 3sin? (3) dK (3.39) 
a 


and this expression is the variance of the log amplitude fluctuations. If 


VAh, with A the EM wavelength, is much greater than unity then the limits 
of integration may be extended from zero to infinity with little error 


[Ref. 18]. Equation 3.39 becomes: 


-8/3 . 2 K7h 
B (0) = mx sin (=*) dK (3.40) 
0 
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and this expression has been evaluated by Tatarski in Chapter 8 of Refer- 
ence 18. In the notation of this calculation the variance of log amplitude 


fluctuations for the layer at a mean altitude of h is given by: 


B (0) = .564 Ah hk”? (C”) (3.41) 


The integral of equation 3.38 is not analytic and a numerical inte- 
gration must be carried out to find the spatial covariance. Normalizing 
equation 3.38 by the variance and carrying out a numerical integration 
yields the function of Figure 5 with p measured in units of VAh. [Ref. 19] 

The intensity spatial covariance is constructed by appropriate shift- 
ing and scaling of the log amplitude covariance as indicated by equation 
3.35. Assuming that the components of the binary source are of equal 


intensity then equation 3.35 becomes: 
C,(p) = 2B (p) + [B (p-eh) + B (p+6h)] (3.42) 
Normalizing each term of this equation separately: 
c(p) = 2b (p) + b.(p+8h) + b (p-6h) (3.43) 


with b,(pt+0h)=B,(pt+0h)/B,(0). A graph of this function for h of ten kilo- 
meters is constructed in Figure 6 with p,h scaled by VAh. 

For experimental implementation the needed covariance functions 
could be computed beforehand and stored for later use. Another alternative 
would be to find a suitable Chebyshev polynomial expansion of the integral 
of equation 3.38. However, the effects of nonideal optics must be considered 
and this implies a numerical calculation of the log amplitude 


autocovariance. 
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Figure 5. The Normalized Log Amplitude Spatial Covariance 
of Fluctuations [Ref. 19] 
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Figure 6. The Normalized Spatial Covariance of Intensity 
Scintillations 
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3. nideal ical Corrections to the A varian 

The actual detector system is coupled to an optical system with the 
real detector array in the exit pupil. Taking the linear systems viewpoint of 
the optical system then it is viewed as a linear transformation of the optical 
input at the aperture. The individual spatial frequencies in the aperture 
are transformed to the exit pupil via the optical transfer function (OTF). 
The OTF of the optical system is defined as the ratio of the output spatial 
power spectrum to the input spectrum, and it is a spatial frequency depen- 
dent complex quantity with the modulus defining the modulation transfer 
function (MTF) and phase, the phase transfer function (PTF). The optical 
system is assumed on axis and since phase shifts in centered optical 
systems occur off axis then the MTF is the quantity of interest. The MTF 
describes the filtering of spatial frequencies by the optical system. For the 
detector array in the exit pupil a spatially modulated input will result in a 
spatially modulated output that is suitably scaled; however, if the MTF falls 
to zero (optical system cutoff) rapidly or is a rapidly varying function of 


spatial frequency then the output spatial spectrum is found from 


F (KO) ie a= MTF) F (K,0) | (3.44) 


tp input 


by the definition of the OTF. The spatial autocovariance then becomes the 
Hankel transform of equation 3.44 as indicated in the previous section. 

The power spectrum of log amplitude fluctuations in the aperture is 
illustrated in Figure 2 and it is noted that this spectrum falls rapidly to zero 
as a function of spatial frequency. Therefore, the assumption of unity MTF 
is quite good for reasonable aperture sizes (greater than approximately 


20cm). The MTF will be assumed constant for the balance of this 
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calculation; however, experimental implementation of this scheme should 
include a measurement of the MTF of the optical system. 

With the assumption of constant MTF over the inertial subrange, 
the scaling between the aperture and exit pupil is derived from geometrical 
considerations. As is apparent from Figure 7 the average intensity in the 
exit pupil is increased by the factor: 


2 
Di / . (3.45) 
Dy 


with Da the aperture diameter, Dr the diameter of the exit pupil and T; the 
overall transmittance of the optical system. For the linear detector array 
illustrated, the scaled array in the aperture is characterized by: 
A,=A,D, fo 3 Pa=PpDp (3.46) 
/ DE Mo 
E 

with A the detector element area, p a linear displacement along the array 
and the subscripts A and E referring to the aperture and exit pupil respec- 
tively. For the remainder of this calculation the optical system coupled to 
the detector is assumed to have unity transmittance and constant MTF over 
the inertial subrange. Under these assumptions the covariances of the 
previous section are assumed valid. 

Given the assumption of an ideal optical system the interaction of 
the electromagnetic field with the detector elements must now be consid- 
ered. The next section discusses the photon counting statistics of the binary 


source profiling method. 
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Scaling Between Aperture and Exit Pupil 
Plane 


40 


IV. PHOTOELECTRON COUNTING STATISTICS 


The intensity scintillations are detected by a linear array of photosensi- 
tive or photoemissive detectors, and the discrete nature of this detection 
process must be considered. The incident electromagnetic flux causes the 
photosensitive surface of a detector element to eject a photoelectron and 
these electrons are then collected and analyzed as the detector element 
signal. Following Saleh [Ref. 25], this derivation treats the electromagnetic 
field classically as a stochastic quantity; and, the photoelectron production 
process is treated semiclassically. This approach does not invoke the 
photon concept and has the advantage of not demanding the full quantum 
electrodynamic treatment of the electromagnetic field. The semiclassical 
approach is compatible with the full quantum analysis of the photoemission 
process provided certain mildly restrictive assumptions hold true. [Ref. 25] 

To proceed with this analysis a certain amount of statistical back- 
ground is required. The definitions and relations used in this calculation 
are summarized in the next section and the notation used is that of 


Reference 25. 


A. STATISTICAL BACKGROUND 

For an experiment that counts photoelectrons the observed data is the 
number of counts (nj,n9,...nm) in the intervals (ty,ty+T) with J equal to 
(1,2,...m). Since the photoemission process is assumed probabilistic, the 


number of counts in disjoint time intervals are assumed independent. 


Therefore, a rate density B(t) is defined: 


41 


(N, (t + At) — N,(t)) 
si (4.1) 
At0 At 


with N7y(t) the number of events in the interval T. This function is also 
known as the arrival rate of a Poisson process [Ref. 26]. 

Given an assumption of independence in disjoint intervals then the 
process described is called a Poisson point process (P.PP) with rate density 
B(t). The number of points (counts) in a single interval (t,t+T) has a Poisson 


density: 


n 


Jeg) =. exp(—W) (4.2) 


with W, the integrated rate of events defined by the relation: 
t +T 


0 
ie | B(t) dt (4.3) 
t 


1] 
And, the time T is the counting time. If B is a constant then the P.PP is 


called homogeneous and equation 4.2 becomes: 





Pu) = St exp (-BT) (4.4) 


n! 


The integer random variable n is also described statistically by the 
moment generating function (mgf). The ordinary moment generating 


function is defined by: 
Q_(s) = (exp Cns)) = 24 _. exp (ns) P(n) (4.5) 


with P(n) given by equation 4.2. Carrying out the indicated summation 


using equation 4.2 gives the following expression for the mgf: 
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Q_(s) = exp | Wexp (-s) -1)] (4.6) 


The usefulness of the mgf is apparent when the term exp(-ns) is expanded 
in a Maclaurin series. The mth term of the expansion evaluated at s =0 is 
the mth ordinary moment of n. Using equation 4.2 and the megf the rela- 
tions between the moments of W and n are easily derived. Several of these 
are listed here for future reference: 


- ordinary moments 


(n)=W a) 
(n2) = W2+W b) (4.7) 
(nyng) = Wi We Cc) 


- central moments 


(6n) = ((n-(n))) = 0 a) 
((8n)*) = W b) (4.8) 
(6n16ng) = 0 Cc) 


The statistics of a non-homogeneous P.PP are completely determined by 
the rate density B(t). However, if B(t) is itself a stochastic function then the 
statistics of the point process depend on the statistics of 8. If B(t) is a 
stationary stochastic function then the point process described by it is called 
a doubly stochastic Poisson point process (DSP.PP). The moments and 
moment generating functions of the DSP.PP are found by averaging the 
corresponding moments and mef's of the P.PP over the different realiza- 
tions of B(t) or the integrated rate W. Thus, the moments of the DSP.PP are 
obtained by averaging equation 4.7 over W. They are listed here: 

- ordinary moments of DSP.PP 
(n) = (W) a) 


(n2) = (W2)+(W) b) 


(4.9) 
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and, conversely: 


(W)={n) a) 
(4.10) 
(W2) = (n2)-{n) b) 
Also: 
(ning) = (WiW92) (4.11) 
- variance of DSP.PP 
((6n)2) = ((6W)2)+(W) (4.12) 


These relations and the results of section III are used to derive the photo- 


electron counting statistics for the binary source method. 


B. THE RATE OF PHOTOEMISSIONS 

A full quantum analysis of the photoemission process reveals the 
following facts [Ref. 25): 

1. The probability of transition from a bound to an unbound state for 
an electron in a photosensitive surface during a time At is proportional to At 
and the instantaneous intensity. 

2. If the intensity is a stochastic function of time then the photoemis- 
sion process is a DSP.PP. 

These conditions hold quite generally if the following two auxiliary 
conditions hold: 

1. The time At, though much greater than the period of oscillation of 
the EM wave, is much shorter than other characteristic time scales of the 
experiment. 

2. The bandwidth of unbound (free) electron states is much greater 


than the EM bandwidth. 
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Under these conditions then the rate density of photoelectric emissions 


can be expressed as: 


bir) = J), a_ dv d°r a (v) I@,t,v) (4.13) 
“Dn 


with v the EM frequency, Ap the detector element area and a(v) the quan- 
tum efficiency at frequency v. 

This expression for the rate density is simplified by making the quasi- 
monochromatic approximation. The quasimonochromatic approximation 
assumes that the center EM frequency is much greater than the EM band- 
width. For optical frequencies this condition is easily satisfied by using an 
appropriate set of colored filters. However, the same thing is accomplished 
by a narrowband detector response. If a more accurate assessment of the 
frequency dependence is required the stellar spectrum is modeled by a 
Planck blackbody spectrum and the explicit frequency dependence of the 
quantum efficiency is used to evaluate the frequency integral of equation 
4.13. However, since an analytic form for the quantum efficiency is not 
readily available this calculation proceeds by invoking the quasimonochro- 


matic approximation and approximating the frequency integral as follows: 
| 2 
Ba,t)=Ava(v)jJ, dr1G,tyv_) (4.14) 
O An O 


with v, the center frequency. 

The spatial integral can similarly be approximated by assuming that 
the area over which the spatial fluctuations are correlated is much larger 
than the detector element area. The rate density then becomes: 

B(r,,t) = Ava(v) Ap Iq.,t,v,) (4.15) 
with r; the geometric center of the detector element. 
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The rate of photoemissions is given by equation 4.3: 
t+T t+T 


W(r,t) = Joceaya = 1 fice dt (4.16) 
t IE 


with y equal to Ava(v,)Ap and letting the frequency dependence be under- 
stood. But, as noted in section III, the intensity is the sum of two indepen- 
dent sources: 

I(r,,t) = I,(x,,t) + 1,@,,t) (4.17) 
Following Fried [Ref. 23], the intensity is expressed as a function of the log 
amplitude as: 


L(x,,t) = F(t) exp (2x,(,,t)) (4.18) 


with I;°(t) the source intensity; and, the exponential term is the atmo- 
spheric modulation. Again, the source term is assumed spatially, though 
not temporally, coherent over the detector element area. Under these 


assumptions the rate is: 
t+T 


Wir,,t) = ¥ Jace exp (2x, (r.,t)) + 1}(texp(2x,(r,t))|dt (4.19) 
t 


The stellar sources are assumed thermal and as is well known thermal 
sources have a coherence time that is the inverse of the optical bandwidth. 
However, the correlation time associated with the atmospheric modulation 
is on the order of one millisecond. This is consistent with the assumption of 
frozen turbulence for times less than or equal to one millisecond [Ref. 20]. 
For any reasonable optical bandwidth the source coherence time is safely 


assumed to be many orders of magnitude less than the correlation time of 
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the atmospheric modulation. Selecting an integration time that 1s much 
less than the atmospheric correlation time and much greater than the 


source coherence time the rate becomes: 


Wir,.t,) = W(t) exp (2x, (r,,t.)) + Wott.) exp(2x,(r,,t,)) (4.20) 


with tj = t+T/2. Note that the atmospheric modulation is considered frozen 
over the integration time. 

Equation 4.20 is used to form the second order moments of the rate and 
these moments are related to the second order moments of the photo- 


electron counts. 


C. COVARIANCE OF PHOTOELECTRON COUNTS 


Forming the normalized spatial covariance of the integrated rate yields: 


(SW(r,,t) 5Wr,,t.)) — (WOr,,t,) W0r,,.)) 
Se ee ee eee ee | (4.21) 


(Wer, ,t.) CWir,,t.))  (W(r,,t.))(W(r,,t,)) 


with the angled brackets denoting an ensemble average. Examining the 
correlation term of equation 4.21 and noting that the atmospheric and 


source modulation are independent then this term becomes: 


(Wor,,t.) Wir,,t.)) = (writ) Xexpl2x,(r t.) + 2x,(r,,t.)]) 


j’ 
+ (Wo(t.) (exp 2x,(r,,t.) + 2x,(r,,t.) 1) + (Wo(t,) Welt.) 


° {(expL 2x, (r,,t.) + 2x,(rort,)|)+(expl 2x, (r,,t,)+2x,(r,.t,)])} (4.22) 


by using equation 4.20. This relation is further simplified by using the 


relation between the intensity and the log amplitude. Following Fried 
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[Ref. 23], the exponential terms of equation 4.22 are related to the log ampli- 


tude covariance by: 
(exp [2x.(r_,t,) + 2x,(r_st)]) = exp 4B, (9) (4.23) 


Since the log amplitude fluctuations are assumed small then this expres- 


sion is expanded in a Maclaurin series to yield: 


exp [4B (p)]~4B_ (p) +1 (4.24) 
1) 1,) 


The various log amplitude covariances implied by equation 4.24 are again 
simplified by appealing to the symmetry of Figure 3 and equations 3.31 to 
3.34. The result is: 


(Wr t.) Wir,,t,)) = [(we(t)) + (Wwett,))] (4B, (p) +1) 
+ 4(WY(t.) W5(t.))[B_(p-@h) + B (p+6h)]} (4.25) 
or, rearranging terms: 
2 
(Wr, tW(r,,t.)) = (Lwoct, #weoctd] ) + [W2(t.)*)+(Wolt,)) ]4B.(p) 
+ 4(W(t.)Wo(t.))[4B (p-Ch) + 4B (p+6h)] (4.26) 


Noting that (exp[2x(r,t)]) equals unity by conservation of energy then the 


normalizing term of equation 4.21 is expressed as: 


(Wr, ,t.))(W(r,,t)) = [(WE(t)) + (wott,))) (4.27) 


by using equation 4.20. Therefore, substituting from equations 4.26 and 4.27 
in equation 4.21 yields the following expression for the spatial covariance of 


the rate: 
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(SW(r, ,t.) SWir,,t.)) ([wi(t,) + wt] ) 
(Wirt, Wert) [ewe «(wee s)]” 
[Cwyet +CWolt)] awit) Ett) 
[(wet))+(wect))] -— [(w°ct)(wet))] 
¢ [B,(p+0h) + B (p-Oh)] - 1 (4.28) 


As stated previously the photoelectron emission process is a DSP.PP for 
a stochastic intensity. The relations between the second order moments of 
W, the integrated rate and n, the number of counts, are given by equations 
4.9 to 4.12. Substituting from these relations in the rate spatial covariance 


of equation 4.28 yields: 
(ante t) Ineryt)) (nit + nytt) — [¢n,(t)) + (ng(t,)] 
(nir,.t)Xnlryt))  [(n,t,))+ (nytt) [nyttp) + (ngtt))T 
[(n,(t,)") o (n,(t)”) - (n,(t.)) a (n,(t.))] 
ee 4B 
[(n,(t)”) + <n, (t)))- 


4(n i(t.)) (n,(t.)) 
[(n,(t,)) + {ngit))) 


This is the normalized spatial covariance of photoelectron counts for a 


xP) 


[B(p-6h) + B (p+6h)] (4.29) 


binary stellar source. 

To proceed further note that a single measurement or counting interval 
T of approximately one millisecond probes a single configuration of the 
atmospheric modulation. This single measurement is clearly insufficient 


to determine the modulation statistics, and since there is only one 


49 


atmosphere to sample then the ensemble averages of equation 4.29 must be 
replaced with a time average. Assuming ergodicity holds, Saleh [Ref. 28] 
points out that long time averages wash out the stochastic nature of the 
intensity fluctuations. To demonstrate this assume that the unmodulated 
source intensity is characteristic of linearly polarized thermal light. Saleh 
derives the probability density of photoelectron counts in terms of the 


parameter N: 


P(n) = ( —_ ) Gs a (4 a (4.30) 


with N the number of coherence times of the source intensity. Apparently, 
a good determination of the atmospheric modulation is made only after 
many coherence times of the source intensity have elapsed. In the limit 


that N goes to infinity equation 4.30 becomes the Poisson distribution: 


we 





Panes 


iY exP (—{n)) (4.31) 


And as is well known [Ref. 26], the following relations hold for the Poisson 


distribution: 


((8n)*) = (n) (4.32) 


Coe Gy aG (4.33) 


Therefore in the limit of long (compared with one millisecond) time aver- 


ages the normalized spatial covariance of counts is: 


00 


2 2 
(Sn(r, t,)n(r4,t;)) (n,) iu (n,) 4 (n,) (n,) 


= 2 4B (p) + 
(n(r,,t.)(n(r,,t.)) [(n, )+(n,)} [(n,) 4 (n.)} 
¢ [B (p-6h) + B (p+6h)] (4.34) 


This expression indicates that the spatial covariance of counts is directly 
dependent on the second order statistics of the atmospheric modulation. 
The source modulation is averaged away due to the extremely short 
(Av 2108Hz) coherence time of the source radiation. 

Since the actual experiment depends on replacing ensemble averages 
with time averages the statistics of the experimentally estimated spatial 
covariance must be explored. The next section examines the counting 


experiment for a single layer Ah. 


D. STATISTICAL ACCURACY 


The spatial covariance of counts: 


= i (4.35) 
(n(r,,t)) (n(ry,t,)) (n(x ,t,)) (n(r,,t,)) 
can be estimated by the expression: 
go = & (4.36) 


A(r,) A(ry) 7 
with; 


2 n(r.,t.) (4.37) 


Z| 


a 1 A 
G = wt rt) n(r,,t_) ; n(r.) = 


ol 


N is the number of counting intervals of length T, and T is expected to be 
much longer than the source coherence time and less than the atmospheric 
modulation correlation time. The estimator g is characterized by its mean 
and variance. The mean or expected value of g is apparently given by equa- 


tion 4.34. 

(n,) +(n,) 4(n,n,) 
—amnneernsps: 15). ((5))) > —_ are 
[(n,)+(n,)) * [n+ dn, )] 


¢ [B (p-@h) + B_(p+6h)] (4.38) 


() = 


However, the variance of g is a difficult quantity to derive. The variance of g 


is defined by the relation: 


var (8) = (8-8) ) (4.39) 


which simplifies to: 


var (g) = = ae eo \] ) (4.40) 


n(r, )n(r,) n(r,)n(r,) 


This quantity is the same as the variance of the estimated spatial correla- 
tion of counts. Saleh evaluates this variance by expanding G and nj about 
their means and retaining only the deviations from these means. With the 
additional assumption that N is very large the variance of g is given by 
equation 7.181 of Reference 25. In the notation of this calculation, the vari- 


ance Of g is: 


2 


var(g) = var(G/ ((n(r,))(n(r,))) +A 


-2(8)(GACr, )) (3° A 
A= 4(g) + fe xa Gees aNotep ate) vr: ve \ (nt, 2 var( ©) | 


(n(r » (n(r,)) 


with the indicated variances the full spatiotemporal variances of the indi- 


(4.41) 


cated quantities. Unfortunately, a complete evaluation of this expression 
requires the time statistics of the atmospheric modulation. The turbulence 
in the atmosphere is typically intermittent in intensity and the statistics of 
intermittency are inaccessible with current theoretical models. However, 
the single scattering assumption allows a lower bound to be placed on the 
variance of g. 

The time statistics of the evolving refractive index inhomogeneities will 
depend on moments of the velocity field of higher order than is represented 
by the structure function (see section II). This implies that the temporal 
variance of the spatioangular correlation of log amplitude fluctuations will 
be a quantity that is fourth order in the log amplitude. By the single 
scattering approximation these moments are approximately zero. There- 
fore neglecting the terms in the variance of g that are of order g2 and 


higher, then equation 4.41 becomes: 


var f 2 9 24g) a 
(g) > 4(g) (ace Mntes) 4(g) [1+ iG Way, 4 (4.42) 


with nj, ng the counts due to source 1 and 2 respectively. The relative error 


is defined as: 


O3 


AZT 
e = 100 [var (g)] is) (4.43) 


OT, 


a [14 .5/[(n,)+(n,)P-] (4.44) 


(() 


A graph of the relative error for (n1) = (ng) =n is given in Figure 8 as a 
function of the spatioangular covariance of counts (g). 

From Figure 8, the relative error diverges rapidly as (3) goes to zero. 
Since observed values of Cn? are very small (10-12 to 10-14cm¥3) [Ref. 12] 
then the relative error induced due to the detection process will dominate 
the counting statistics. Clearly some error reduction scheme is required to 


extract the spatial covariance of counts. 


o4 


2x10 





Figure 8. Relative Error Versus Spatioangular Covariance 
of Counts for Various Values of Incident Intensity 
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V. LEAST SQUARES ANALYSIS 


The analysis of the preceding section indicated the dominance of the 
detection uncertainty in the determination of the spatioangular covariance. 
This large error is ameliorated by using the redundant information 
provided by the multiple detector elements of the array (Figure 3). Least 
squares analysis is employed to use this redundant information to reduce 
noise. 

Most of the analysis to this point has applied to a single layer at a mean 
altitude h. However, equation 4.38 is easily generalized to any layer with a 


mean altitude of h;: 


(oy) (og) 


C (p,,h,) = 4B (p.h)4 
©) En, )4(n,)] 
4(n, Xn,) | , 
— . 7 x42 LB. (p.-6h,) + B(p,+6h.) (5.1) 
[(n, )+4n,)] 


Note that the log amplitude covariances are now calculated at the centers of 


the detector elements pj. With the previous assumption of a small detector 
element size this is a reasonable approximation. Rewniting the spatio- 
angular weighting functions to reflect the explicit dependence on the 


refractive index structure parameter yields: 
2 1 
B, (ph) = (C*(h,)) B (ph) (5.2) 


and equation 4.38 is generalized to: 
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(n,)'+n,) 
[(n,)+(n,)] 


Z 1 
C,(p,,h,) = (C*(h.)) { 4B\(p,,h;) + 


4{n,)(n,) 


7. [B,(p;-6h,) + B(p,+8h;)} } 


Z 
[(n, +(n,)] 
= (C;(h,)) R(p,,h,) (5.3) 


Using this equation, and assuming the fluctuations induced by non-over- 
lapping layers are uncorrelated, then spatioangular covariance of counts 


observed in the aperture due to multiple layers is given by the summation: 


2 
CR(p,) = 2 C,(p,h,) = 2 (c?(h,)) Rip,.h,) (5.4) 


1 


This summation is equivalent to the matrix equation: 


[cr,] = [R,,] [(c2)] (655) 
with CR, the experimentally measured spatial covariance of counts, Cy?; 
the undetermined refractive index structure parameter profile and Rj, the 
weighting functions described above. The undetermined Cy” are found by 
minimizing the quantity: 


| [CR, ] - [R,,] [(c,),] | = minimum (5.6) 


as a function of C,*. The solution to this is well known and is given by: 


[o? J=((R, J TR, [R,) [eR] 6.7) 


with the superscript * denoting the Hermitian conjugate. This process is 


formally known as least squares analysis and the C,? vector obtained is the 


o7 


least squares approximation of the structure parameter profile. A more 
explicit expression for the profile is obtained by considering the experimen- 
tally derived spatioangular covariance. 


The experimental covariance is given by the matrix: 


[CR] = [a(r,)n@,) =i) (5.8) 


Assuming the angular separation of the binary source is small then the 
expected values of these matrix elements are a function only of the separa- 
tion between detector elements. Therefore, elements along diagonals 
parallel to the main diagonal are expected to be equal in the absence of 
noise. Using this fact, the spatioangular covariance of counts as a function 
of detector element separation is formed by summing along these diagonals 


as indicated here. 


n(r,)n@,) nj n(r5)-.- n@ na 
nUGeyNoey). > NGaye S Se 


N(rg)A(ry) > Alrg)Arg) +, 


ae SS (5.9) 


A(r,)ACr,) * A(x yo )A(rq9) 


[CRj.] > [CRy] by 





Since the number of terms along these diagonals are not equal, then the 
elements of this sum must be normalized by the number of terms to 
preserve the relationship of equation 5.4. Therefore, the elements of the 


measured covariance vector are written: 
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—1 
CR, =(2(10-In-m!)) 2. [CR +CR__] (5.10) 
n=1,10-k 
m=n+k 
and, the structure parameter profile is determined by solving the matrix 


equation 
[2 (R, R,,)] [e?m)] = [2 (R, , CR,)] (5.11) 


Which results from applying least squares analysis with the parameters 


Cn2(hj) the undetermined coefficients. The profile is now determined since 


the detector element separation and altitude are connected by the relation 


Sul., 


og 


VI. CONCLUSIONS 


This thesis has explored the theoretical basis for profiling the refractive 
index structure parameter using the spatial covariance of binary source 
intensity scintillations. This analysis derived the power spectrum and 
associated spatial covariance of intensity scintillations caused by atmo- 
spheric refractive turbulence. The expected spatial covariance of 
photoelectron (photon) counts and counting statistics were derived using 
the intensity covariance functions. 

No rigorous analysis of the binary source technique as implemented by 
several French teams [Refs. 12-15] is currently available in the open litera- 
ture; and, the analysis of this thesis reveals several interesting features not 
previously noted. 

The power spectrum and spatial covariance used by the French teams 
are essentially correct; however, the error analysis of the binary source 
technique is incomplete. A consideration of the photoelectron counting 
statistics indicates that very large relative errors of 200 to 2,000 percent are 
expected in any determination of the spatial covariance function. This 
large error is not significantly improved by increasing the counting rate 
above unity as indicated by Figure 8. Previous experimental work [Refs. 
12-15] was limited to bright (magnitude 2) stellar binaries. However, an 
elementary calculation of the available photons per counting time 
(approximately one msec) coupled with the relative insensitivity of the 


binary source technique to increased counting rates indicates that stellar 
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binaries of magnitudes three to four are also suitable for detector apertures 
of 30 to 60 centimeters. 

The large relative error expected in an experimental determination of 
the spatial covariance necessitates the use of a smoothing algorithm and 
the least squares technique was selected. The advantage of least squares 
analysis is that it accomplishes the two-fold purpose of data smoothing and 
refractive index structure parameter profiling. The basis selected for the 
least squares analysis is the nonorthogonal set of theoretically derived 
spatial covariance functions for a series of atmospheric layers of increasing 
attitude. The weighting coefficients derived using the least squares algo- 
rithm are found to be the required refractive index structure parameters for 
the atmospheric layers. 

The expected improvement in the relative error remains an open-ended 
question. The error analysis depends on the intermittency of atmospheric 
turbulence and as indicated in Section IV, the error analysis in this thesis 
constitutes a lower bound for the expected error. Clearly, further work on 
the intermittent aspects of atmospheric turbulence is indicated. 

Despite the large errors expected in the determination of the spatial 
covariance function the binary source technique has several clear 
advantages. As mentioned in the Introduction, the altitude resolution of 
the binary source technique is superior to that of existing methods. The 
implementation of the least squares algorithm is straightforward and 
results in a relatively unambiguous determination of the structure 
parameter profile. Modest aperture sizes (30 to 60 cm) are adequate for the 


detector system. The greatest disadvantage of the technique is the paucity 


of suitable stellar binaries. However, by using stellar sources of up to 
magnitude four this drawback becomes minor. 

Based on this analysis, experimental implementation of this technique 
should proceed, and the experimental profiles obtained should be compared 
with insitu measurements or integrated profiles such as that produced by 


stellar scintillometers [Ref. 16]. 
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